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Abstract 

Approximating ground and a fixed number of excited state energies, or equivalently low order Hamil¬ 
tonian eigenvalues, is an important but computationally hard problem. Typically, the cost of classical 
deterministic algorithms grows exponentially with the number of degrees of freedom. Under general 
conditions, and using a perturbation approach, we provide a quantum algorithm that produces estimates 
of a constant number j of different low order eigenvalues. The algorithm relies on a set of trial eigenvec¬ 
tors, whose construction depends on the particular Hamiltonian properties. We illustrate our results by 
considering a special case of the time-independent Schrodinger equation with d degrees of freedom. Our 
algorithm computes estimates of a constant number j of different low order eigenvalues with error 0(e ) 
and success probability at least |, with cost polynomial in I and d. This extends our earlier results on 
algorithms for estimating the ground state energy. The technique we present is sufficiently general to 
apply to problems beyond the application studied in this paper. 

The final publication is available at Springer via http://dx.doi.org/10.1007/slll28-015-0927-y 


1 Introduction 

Computing eigenvalues of Hamiltonians with a large number of degrees of freedom is a very challenging 
problem in computational science and engineering. Hamiltonian eigenvalues give the system energy levels, 
corresponding to the ground and excited states. For example, one of the most important tasks in chemistry 
is to calculate the energy levels of molecules, which is required for predicting reaction rates and electronic 
structure, and which, in particular, depends principally on the low order energy levels. The best classical 
algorithms known for such problems have cost that grows exponentially in the number of degrees of freedom 
[22] . Therefore, efficient quantum algorithms would be an extremely powerful tool for new science and 
technology. 

On the other hand, there are a number of recent results in discrete complexity theory suggesting that 
many eigenvalue problems are very hard even for quantum computers because they are QMA-complete m 
15511551 15], However, discrete complexity theory deals with the worst case over large classes of Hamiltonians. 
It does not provide methods or necessary conditions determining when an eigenvalue problem is hard. In 
fact, there is a dichotomy between theory and practice. As stated in [55], “complexity theoretic proofs of 
the advantage of many widely used classical algorithms are few and far between.” Therefore, it is important 
to develop new quantum algorithms and to use them for solving eigenvalue problems for which quantum 
computing can be shown to have a significant advantage over classical computing. 

In [29] we developed an algorithm and proved strong exponential quantum speedup for approximating 
the ground state energy (i.e., the smallest eigenvaue) of the time-independent Schrodinger equation under 
certain assumptions. In [30] we explain why this problem is different from the QMA-complete problems of 
discrete complexity theory. In |28j we relaxed an important assumption of |29j and extended our results 
to the ground state energy approximation for the time-independent Schrodinger equation with a convex 
potential. 
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An important advance would be to obtain analogous results for approximating excited state energies 
under weakened assumptions. The techniques we have used previously for the ground state energy do 
not extend to excited state energies. Similarly, in computational chemistry, for instance, Holienberg-Kohn 
density functional theory (DFT) is strictly limited to ground states [141 T8] . There are other flavors of DFT 
that may provide approximations of excited state energies. However, in general, approximate methods in 
computational chemistry often succeed in predicting chemical properties yet their level of accuracy varies 
with the nature of the species and may fail in important instances; see mm and the references therein. 
Obtaining conditions allowing one to approximate excited state energies with a guaranteed accuracy and a 
reasonable cost would provide a valuable insight into the compexity of these problems. 

In this paper we present an entirely new approach for approximating a constant number of low order 
eigenvalues. At the same time we relax some of the assumptions of our previous work for approximating 
the ground state energy [29] [28]. We will discuss these papers in Section l2Tl Using the properties of 
our eigenvalue problem, we construct a set S of trial eigenvectors. This set contains vectors that overlap 
sufficiently with the unknown eigenvectors corresponding to the eigenvalues of interest. Then, these vectors 
can be used as initial states in quantum phase estimation (QPE) to produce eigenvalue estimates with a 
reasonably high (i.e., not exponentially small) success probability. The elements of S are known eigenvectors 
of a slightly perturbed problem. It is important to select the perturbation carefully so that the elements of 
S can be prepared efficiently on a quantum computer. In principle it is difficult to determine exactly which 
eigenvectors of the perturbed problem sufficiently overlap with the unknown eigenvectors of interest. Thus, 
our construction of S generally contains more elements than are absolutely necessary. Our algorithm runs 
QPE repeatedly with each element of S as initial state. We show that carefully selecting a constant number of 
the smallest measurement outcomes leads to estimates of the desired eigenvalues with a reasonable probability 
and cost, as long as the size of S is not exponentially large in the problem parameters. By reasonable cost 
we mean that the algorithm uses a number of qubits and quantum operations which is polynomial in the 
problem parameters. By reasonable success probability, we mean a probability p that is bounded from below 
by a constant, e.g. p > |. Unless the success probability of an algorithm is exponentially small in the 
problem parameters, it can be boosted to become arbitrarily close to 1 using a number of repetitions that 
is also polynomial. We remark that the selection of the perturbation of the Hamiltonian impacts the size of 
<S, and hence the cost of our algorithm, and is an important consideration. 

We illustrate our results by considering the time-independent Schrodinger equation under weaker as¬ 
sumptions than those of [29] [28], as we explain below. For this problem, the cardinality of the set S of 
trial eigenvectors turns out to be polynomial in the number of degrees of freedom d. We derive cost and 
probability estimates for approximating a constant number of low order eigenvalues. Indeed, for accuracy 
0(e) the cost and the number of qubits of our algorithm is polynomial in d and A More precisely, we 
consider the eigenvalue problem 

(^-Ia + V^*(x) = E*(x) x g Id = (0, l) d , (1) 

T(z) = o xedi d , (2) 

where A denotes the Laplacian and $ is a normalized eigenfunction. Here all masses and the normalized 
Planck constant are set to one, and we assume the potential V is a smooth and uniformly bounded function 
as we will explain later. 

Our problem is to compute a constant number j of estimates 

E 0 < Ei < ... < Ej- 1, 

each approximating a different low order eigenvalue with error 0(e) and high probability. Eq is the estimate 
of the ground state energy (i.e. the smallest eigenvalue). In general, the eigenvalues may be degenerate 
with unknown multiplicities. Moreover, eigenvalues may be clustered in balls of radius 0(e). We call such 
eigenvalues e- degenerate. It is reasonable to assume that it is not necessary to produce estimates for every 
single e-degenerate eigenvalue, and some of them can be omitted. 

Such eigenvalue problems can be solved by suitably discretizing the continuous operator (Hamiltonian) 
to obtain a symmetric matrix whose low order eigenvalues approximate those of the continuous problem, 
and then by approximating the matrix eigenvalues. Eigenvalue problems involving symmetric matrices are 
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conceptually easy and methods such as the bisection method can be used to solve them with cost proportional 
to the matrix size, modulo polylog factors HU. The difficulty is that the discretization leads to a matrix of 
size that is exponential in d. Hence, the cost for approximating the matrix eigenvalue is prohibitive when 
d is large. In fact, a stronger result is known, namely the cost of any deterministic classical algorithm 
approximating the ground state energy must be at least exponential in d, i.e., the problem suffers from the 
curse of dimensionality |27] . It is important to point out that different approaches may lead to different matrix 
eigenvalue problems that have varying degrees of difficulty. For instance, in quantum chemistry, the first and 
second quantization approaches for computing energies of the electronic Hamiltonian, as described in m , 
lead to completely different matrices with different notions of degrees of freedom. Moreover, discretizations 
of certain problems in physics may lead to eigenvalue problems for stoquastic matrices, that some believe 
are computationally easier to solve [7]- 

It is worth noting that in certain cases quantum algorithms may be able to break the curse of dimen¬ 
sionality by computing e-accurate eigenvalue estimates with cost polynomial in e _1 and d. This was shown 
in .29,121] where we saw that for smooth nonnegative potentials that are uniformly bounded by a relatively 
small constant, or are convex, there exists a quantum algorithm approximating the ground state energy with 
relative error 0(e) and cost polynomial in d and e” 1 . 

It is important to investigate conditions for the potential V beyond those of [2UI2S1 HO [3D]. In this 
paper we pursue this direction. As we indicated, we give a general algorithm for low order eigenvalues, and 
then apply it to the time-independent Schrodinger equation where V is smooth and uniformly bounded by 
a constant. The algorithm has cost polynomial in e _1 and d , regardless of the size of the bound. We exhibit 
the resulting quantum algorithm, its cost, and success probability. The technique that we have developed 
can be applied to other eigenvalue problems as well. 

We summarize the contents of this paper. In Section 2 we define our eigenvalue problem. We also review 
classical and quantum algorithms for eigenvalue problems. We discuss the limitations of classical algorithms, 
and how they may be overcome by quantum algorithms. We specify rather general conditions and provide 
a quantum algorithm which computes a constant number of approximations to low order eigenvalues, i.e. 
low order excited state energies (including the ground state). We explain how to construct a set S of trial 
eigenvectors for our algorithm using a perturbation approach. In Section 3, we study the overlaps between 
the trial eigenvectors and the unknown eigenvectors corresponding to the eigenvalues of interest. We provide 
lower bounds for the overlaps, and show how they depend on the cardinality of S. In Section 4 we illustrate 
our results by considering a special case of the time-independent Schrodinger equation. This allows us to 
present specific estimates for the cost of our algorithm and its success probability, which we state explicitly 
in Theorem [I] Finally, we summarize our results in Section 5. 


2 Problem Definition 

In this section, we introduce the eigenvalue problem in its most general form to emphasize that our approach 
applies under very broad conditions. In later sections, we will make more assumptions in order to show 
specific results. 

We consider an eigenvalue problem for a self-adjoint operator L with a discrete spectrum. We will provide 
more details about L below. Let 

-®(o) < -®(i) < ) < ... (3) 

be its eigenvalues ignoring multiplicities, which we call the energy levels of L. Suppose we want to estimate 
the lower part of the spectrum with accuracy 0(e). Since any two distinct eigenvalues of L can be arbitrarily 
close to each other, any algorithm that approximates the lower part of the spectrum with accuracy e cannot 
be expected to distinguish between all eigenvalues E^ ^ Em with | E^ — E^\ = O(e). We have called 
such eigenvalues e- degenerate. So the goal is to obtain an algorithm whose output will be j numbers 

Eq < Ei < ... < Ej-i (4) 


satisfying with high probability the following conditions: 

Cl For every i ^ k G {0, ... ,j— 1}, there exist E ( Si ) ^ ^(s fe ) such that \E^ Si ^—Ei\ = O(e) and |£!( Sfc ) — Ek\ = 
O(e), i.e. different outputs are approximations of different eigenvalues with error O(e), respectively. 
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C2 If \E i+ i — Ei\ = ta(e), there is no eigenvalue E of L satisfying Ei < E < E i+ i and min(|£^ i+ i — 
E |, | Ei — E\) = w(e)LU Thus the algorithm doesn’t miss (or skip) any eigenvalues in the lower part of 
the spectrum unless they are O(e) apart, i.e., e-degenerate. 

Clearly, if e is sufficiently small such that the eigenvalues of L are well-separated, then the algorithm produces 
approximations with error O(e) of the j smallest distinct eigenvalues. 

Assume that L° is a self-adjoint operator defined on a separable Hilbert space, V is a symmetric operator 
whose domain contains the domain of L°, and that L = L° + V is self-adjoint on the domain, D(L°), of L°. 
We also assume that L° and L have discrete spectra and that the eigenspaces associated with each eigenvalue 
are finite dimensional; see e.g. pm DU E2]. In the general case, selecting the partition of L to L° and V is 
not trivial and may significantly affect the problem complexity; we do not deal with this problem here. Our 
discussion in this section applies equally well to Hermitian matrices. 

Let 

a < Eq < Ei < ... < Ei < ... ( 5 ) 

be the eigenvalues of L indexed in non-decreasing order, where a is a given lower bound. Ignoring possible 
eigenvalue multiplicities we have a strictly increasing subsequence of eigenvalues which we denote by 


E(o) < E(i) < ... < E(.i) < ... ( 6 ) 

Similarly we denote by 

E§<E?< ... < E° < ... (7) 

the eigenvalues of L° indexed in non-decreasing order, and by 

E° 0) < E° {1) < ... < E° {i) < ..., (8) 

the eigenvalues of L° ignoring multiplicities. Assume we know all the eigenvalues and eigenvectors of L°. 
Often this is a reasonable assumption. For example, this is true for the eigenvalues and eigenvectors of the 
Laplacian L° = — A defined on the d-dimensional unit cube with Dirichlet or Neumann boundary conditions. 

We wish to estimate the low order eigenvalues of L. By low order we mean that j in equation fif is a 
constant. Intuitively, we expect a “small” and suitably well-behaved perturbation to have a proportionately 
“small” effect on the eigenvectors and eigenvalues of L°. Algorithms solving this problem can take advantage 
of the known eigenvalues and eigenvectors of L°. 

2.1 Background: Classical and Quantum Algorithms 

We briefly review algorithms for eigenvalue problems. Recall that we are interested in problems for which 
quantum computing can be shown to have a significant advantage over classical computing with performance 
guarantees in terms of accuracy and speed. Hence, we do not consider empirical approaches or heuristic 
eigenvalue algorithms. 

Algorithms approximating eigenvalues use a discretization of L to obtain a matrix eigenvalue problem. 
For example, when L is a differential operator, one can use a finite difference discretization [23], or a finite 
element discretization El E]- In particular, for the time-independent Schrodinger equation specified by 
equations © and ©, a finite difference discretization has been used in [22; [23]. Since L is self-adjoint, the 
resulting matrix is symmetric. 

Matrix eigenvalue problems have been extensively studied in numerical linear algebra, and there are 
classical algorithms for approximating one, or some, or even all of the eigenvalues and/or the corresponding 
eigenvectors of a matrix nn ©a ei m- Examples of such algorithms include the power method, inverse 
iteration, the QR algorithm, and the bisection method for symmetric matrices. In particular, the bisection 
method for symmetric matrices can compute the eigenvalues that lie within a given range. In general, the 
known eigenvalues of L° can be helpful in computing such a range for the j eigenvalues of L that we consider 
in this paper. Typically, a symmetric matrix will be transformed to Hessenberg form, which is a tridiagonal 
matrix [11 4 Sec. 4.4.7]. Then the bisection method is applied to the latter matrix [TT] Sec. 5.3.4]. The 

1 For functions /,g > 0 defined on R_|_, the notation /(e) = ui(g(s)) means that for any M > 0, arbitrarily large, we have 
/(e) > Mg(e) for sufficiently small e. 
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cost of this procedure, even if the original matrix is dense, is a low degree polynomial in the matrix size 
and log -, where e is the desired accuracy. Therefore, for matrices of moderate size, approximating the low 
order eigenvalues can be done at a reasonable cost. On the other hand, the costs of the above algorithms 
are bounded from below by a quantity that is at least proportional to the matrix size, even if the original 
matrix is sparse. Hence, the eigenvalue estimation problem becomes hard when the matrix size is huge. 

Observe that the discretization of the operator L must be sufficiently fine so that the eigenvalues of L 
which are of interest are approximated by eigenvalues of the resulting matrix within the specified accuracy 
e. This increases the matrix size. For example, in the estimation of the ground state energy (smallest 
eigenvalue) of the time-independent Schrodinger equation, equations © and d2j), with V uniformly bounded 
by 1, the finite difference discretization on a grid will yield a matrix of size m d x m d , m = log 2 — 1 (29] , 
This means that the cost of the matrix eigenvalue algorithms mentioned above is bounded from below by a 
quantity proportional to (-) , i.e. the cost grows exponentially in d. In [29] , the Laplacian was discretized 
using a 2d + 1 stencil on a grid, and V was discretized by evaluating it at the grid points. 

To approximate the ground state energy of the problem specified by equations © and © in the worst 
case with (relative) error e, assuming V and its first-order partial derivatives are uniformly bounded by 1, 
and the function evaluations of V are supplied by an oracle, a much stronger result holds. The complexity 
(i.e., the minimum cost of any classical algorithm, and not just the eigenvalue algorithms mentioned above) 
is bounded from below by a quantity proportional to e~ d as ds —>■ 0[27j[29]. So unless d is moderate, the 
problem is very hard and suffers from the curse of dimensionality. In [30], we elaborate on this lower bound. 
The same complexity lower bound applies to the approximation of low order eigenvalues under the same, 
or more general, conditions on V. Finally, we point out that the complexity of this problem in the classical 
randomized case is an open question. 

We now turn to quantum algorithms. There is a well-studied quantum algorithm, quantum phase esti¬ 
mation (QPE) [I][26]l4], which can be used to approximate eigenvalues of a Hamiltonian H. More precisely, 
the algorithm approximates the phase corresponding to an eigenvalue of a unitary matrix, which in our case 
is e~ lH . QPE is efficient if two conditions are met. The first condition is that simulating a system evolving 
with Hamiltonian H can be done efficiently, i.e. we can approximate ( £ 1, accurately with low 

cost. The second condition requires that we are given a relatively good approximation of an eigenvector 
corresponding to the eigenvalue of interest. In addition, one should be able to implement this approximation 
as a quantum state efficiently. The approximate eigenvector is used to form the initial state of QPE. We 
also remark that QPE uses the (quantum) Fourier transform as a module. The Fourier transform can be 
implemented efficiently on a quantum computer [261 . 

We discuss the two required conditions for QPE further. Simulating the evolution of a system under a 
Hamiltonian H appears to be a difficult problem for classical computers when the size of H is large. As 
proposed by Feynman |12j . quantum computers are able to carry out such simulation more efficiently in 
certain cases. For example, Lloyd [21] showed that local Hamiltonians can be simulated efficiently on a 
quantum computer. About the same time, Zalka |43[ 142] showed that many-particle systems can be also be 
simulated efficiently on a quantum computer. Later, Aharonov and Ta-Shma [2j generalized Lloyd’s results 
to sparse Hamiltonians. Berry et. al. [6j extended the cost estimates of [2- The results of |6j were in turn 
improved by Papageorgiou and Zhang in [31] . Although there has been more work on quantum Hamiltonian 
simulation since then, the approach of HIS] suffices for our discussion. These papers assume that H is given 
by a black-box (or oracle), and that H can be decomposed efficiently by a quantum algorithm, using oracle 
calls, into a finite sum of Hamiltonians that individually can be simulated efficiently. In this paper, where 
L — L° + V, we assume that the Hamiltonians resulting from the discretizations of L° and V can be simulated 
efficiently on a quantum computer. Their sum, i.e. the Hamiltonian obtained from the discretization of L, 
can be simulated efficiently using splitting formulas such as the Trotter formula, the Strang splitting formula, 
or Suzuki’s high-order splitting formulas. Simulation cost estimates are shown in m- 

Moreover, since we know the eigenvalues and eigenvectors of L°, in certain cases one might be able to 
simulate its evolution explicitly without relying on an oracle. For instance when L° = —A, a quantum 
algorithm and circuit implementing the evolution of the discretized Laplacian, is shown in [8]. That paper 
deals with the solution of the Poisson equation with Dirichlet boundary conditions. The efficient simulation 
of the discretized Laplacian is achieved by diagonalizing it using the quantum Fourier transform. The efficient 
simulation of the discretization of V (which is a diagonal matrix) is achieved using an oracle and quantum 
parallelism. 
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We now turn to the second requirement of QPE, namely the availability of a good approximate eigen¬ 
vector. QPE will produce an estimate of the eigenvalue A (or more precisely, an estimate of the phase 
(f> £ [0,1) corresponding to A through A = e 2 "^) with success probability proportional to the quality of the 
approximate eigenvector mm- If the eigenvector providing the initial state of QPE is known exactly, the 
parameters of QPE can be set so its success probability is arbitrarily close to 1 [2B]. If, on the other hand, 
we use an approximate eigenvector, the success probability is reduced proportionally to the square of the 
magnitude of the projection of the approximate eigenvector onto the actual eigenvector (i.e. the square of 
the overlap between the two vectors) [T]. As long as this overlap is not exponentially small, QPE is efficient. 

For example, in [28) we show quantum algorithms that meet the two requirements of QPE and 
approximate the ground state energy for special cases of the time-independent Schrodinger equation, as 
specified in equations © and ([2j). In j2f)] . we assumed V and its first-order partial derivatives are uniformly 
bounded by 1. In that paper we show a quantum algorithm estimating the ground state energy with cost 
polynomial in - and d. In [30], we explain why quantum algorithms have a significant advantage over classical 
algorithms solving this problem in the worst case. In the second paper [28j . we extend the results to a different 
class of potential functions V, namely convex functions uniformly bounded by an arbitrary constant C > 1 
with first-order partial derivatives uniformly bounded by a constant C'. Under these assumptions, we derive 
a multistage quantum algorithm for estimating the ground state energy. The convexity of U, along with a 
recent result 0] concerning the fundamental gap (i.e., the difference between the first two eigenvalues) of 
Schrodinger operators, allows us to use a number of stages that is polynomial in d, with each stage having 
cost polynomial in - and d , so that the overall algorithm is efficient. Therefore, for special cases of the 
time-independent Schrodinger equation, quantum algorithms vanquish the curse of dimensionality. 

We remark that obtaining a good approximate eigenvector required for QPE is a particularly difficult 
task, in general, when the matrix size is huge. Things are complicated further if one needs a number of 
different approximate eigenvectors, in order to use QPE to approximate the j > 1 lower order eigenvalues. 
We overcome this difficulty for L = L°+V using the known eigenvalues and eigenvectors of L°, and properties 
of V, as we discuss below. 

2.2 Algorithm 

2.2.1 Algorithm: Idea 

Our goal is to use QPE to estimate j low order energy levels of L. For this, we need relatively good 
approximations of the corresponding eigenvectors. Since L = L° + V (i.e. L and L° differ by the perturbation 
V ), and since we know the eigenvalues and the eigenvectors of L°, we can use them to obtain the necessary 
approximate eigenvectors. We indicate how this can be done. For simplicity and notational convenience, 
we do not distinguish between operators and their matrix discretizations in the rest of this section, since 
it is not important for the moment. Let E denote one of the low order eigenvalues of L (see, equation 
©) that we wish to estimate. Intuitively, we expect a “small” and suitably well-behaved perturbation to 
have a proportionately “small” effect on the eigenvectors and eigenvalues of L°. Let u be an arbitrary 
unit vector belonging to the eigenspace associated with E. Then there exists an eigenvector u° k of L°, 
similarly corresponding to a low order eigenvalue, that has an overlap (magnitude of projection) with u that 
is non-trivial, i.e. | | is not extremely small, as we will see later. 

A simple illustration of this idea is to imagine an L° with a symmetric ground state, and a perturbation 
V that is relatively asymmetric. In such a case, the ground state of L may overlap primarily not with 
the ground state eigenvector of L° but with an excited state. For example, the ground state of the one¬ 
dimensional harmonic oscillator is a (spatially) even function about the center of the well, and successive 
eigenstates are alternately odd and even. An odd perturbation of sufficient size will cause the ground state 
to become approximately odd, and project increasingly onto an odd unperturbed eigenstate. 

An important idea in this paper is to form a collection S of the eigenvectors of L° that correspond to 
eigenvalues of L° that satisfy a certain property, which we specify in the next subsection. The goal is to have 
at least one element in S that has a reasonable overlap with a vector in the eigenspace corresponding to Eu\ , 
for each i = 0,1, ...,j — 1. We call S the set of trial eigenvectors. We will use each one of the elements of 
S repetitively as initial state in QPE, running QPE multiple times, to obtain a sequence of approximations 
that will lead us to estimates of each E^. 
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Let us briefly discuss the idea for constructing S. At one extreme, one could take S to be all of the 
eigenvectors of L°, because not all of them have a negligible overlap with the eigenvectors of L corresponding 
to the eigenvalues of interest. However, then the size of S can be huge. To limit |<S|, we select eigenvectors 
of L° that correspond to eigenvalues that do not exceed a certain bound. Roughly speaking, we will be 
excluding eigenvalues of X° that correspond to energies grossly exceeding the energies of L that we wish to 
estimate. This idea is made precise in equation (O in next section. 

The cardinality of S depends on the eigenvalue distribution of L°. If the cardinality of S is not pro¬ 
hibitively large, and if we can discretize its elements and efficiently prepare the corresponding quantum 
states, then we can run QPE repeatedly for the all elements of S to produce an estimate of among its 
different outputs with a sufficiently high probability, i = 0,1,..., j — 1. This probability can be boosted to 
become arbitrarily close to 1 using further repetitions of the procedure. We remark that the cardinality of S 
depends on the distribution of eigenvalues of L° and the properties of V. Observe that detecting the desired 
estimates Eq, ..., Ej-i from the outcomes obtained from the different runs of QPE is not a trivial task, and 
we will show how this is accomplished. 

2.2.2 Algorithm: Description 

Let V be such that ||R||lo := sup{||Vw|| : u £ D(L°), ||u|| = 1} < oo uniformly in d. Assume we are given 
(or we have derived) c > 1 and B a sufficiently large upper bound on the lower part of the spectrum of L 
which is of interest o Consider the set of indices 

l:={i:E°-B>c\\V\\ L o}^«). (9) 

We define S to be the set of eigenvectors of L° that correspond to eigenvalues E® with i ^ X; in the case 
of degeneracy, it suffices to select any basis of the degenerate subspace. By constructing S in this way, we 
are guaranteed that at least one of its elements will overlap sufficiently with an element of the degenerate 
eigenspace corresponding to each Eu\, for * = 0,1, ...,j — 1. We will show that the magnitude of this overlap 
is bounded from below by a positive constant. 

The following is an overview of our quantum algorithm for approximating j = 0(1) low order eigenvalues 
of the operator L = L° + V. Algorithm 1 deals with the special case of approximating the ground state energy 
2?(o) • This algorithm illustrates our idea of using a set of trial eigenvectors to approximate an eigenvalue of 
L. Algorithm 2 computes the sequence of approximations E±, E^, ..., £y_i, where each Ei is computed using 
the values Eq through £)_i. Thus the overall procedure consists of iterating Algorithm 2 until we obtain 
the j desired estimates of equation a. 

Let us pretend for the moment that L° and L are N x N matrices; we do this for notational convenience. 
In later sections, when we consider specific instances of L° and L, we will show how to discretize them and 
obtain symmetric matrices such that each of the low order eigenvalues of these matrices approximates the 
corresponding eigenvalue of the respective operator with error proportional to s. 

Our algorithms are based on QPE and require two quantum registers. The first register ( top register) 
contains sufficiently many qubits t to guarantee the required accuracy 0(e ) in the results with a reasonable 
success probability for QPE. The second register ( bottom register) contains the necessary number of qubits 
to hold an approximate eigenstate. 

Algorithm 1. Description - Ground State Energy: 

1. Define S , the set of trial eigenvectors, to be all eigenvectors of L° that correspond to eigenvalues E® 
i X as defined in equation ©■ We denote these eigenvectors by vfj, for k = 0,1,.., |«S| — 1. 

2. Set k = 0. 

3. Prepare the initial quantum state |0)® 4 The value of t is chosen so that QPE, with relatively high 
probability, produces outcomes leading to energy estimates with error O(e). 

4. Perform QPE with initial state |O) 0t |u°) using the unitary matrix U = e lA ^ R . A = L if L non-negative 
definite, and otherwise A = L — al, where a is a lower bound to the minimum eigenvalue of L as we 

2 We give an explicit construction for B in equation CCD. 
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have assumed in the previous section. The parameter a is assumed to be known; see equation ©. 
Nevertheless, even if a is not known, it is often possible to obtain a convenient estimate of a using the 
eigenvalues of L° and the properties of V. The goal is that A is a non-negative definite matrix. The 
parameter R is an upper bound to the spectral norm of A , which can be obtained using the eigenvalues 
of L° and the properties of V. The purpose of R is to ensure that the resulting phases will lie in the 
interval [0,1). 

5. Measure the first t qubits, which give the result of QPE, and store the resulting value classically. We 
assume that the measurement outcomes are truncated to b bits and we obtain non-negative integers in 
the range {0,..., 2 b — 1}, where b < t. The role of the extra qubits to = b — t is to increase the success 
probability of QPE. 

6 . k <— k + 1 . 

7. Repeat steps 3-6 while k < |<S|. 

8. Repeat steps 2-7 r many times, where r is a number precomputed to ensure with high probability that 
the stored results after r runs contain an estimate of Eq. The value of r depends on the problem at 
hand. In Section [4] we derive r for a particular application. 

9. Take the minimum value of the stored measurement outcomes, mark it as selected, and convert it to 
an eigenvalue estimate Eo of Eq using the values of a and R in the definition of A. 

10. Output Eq- 


Note, the purpose of step 7 is to run QPE |<S| many times, once with each |u°) £ S as input, because we do 
not know which of the elements of S has the largest overlap with the unknown ground state eigenvector, and 
the success probability of each run depends on this overlap. Since the largest overlap between the elements 
of S and the unknown eigenvector may not be sufficiently large so that the resulting success probability of 
the algorithm is bounded from below be a constant, say |, the purpose of step 8 is to repeat the entire 
procedure r many times to boost the success probability of computing Eo correctly. 

The following iterative algorithm extends Algorithm 1 to compute the sequence of approximations 
Ei,..,Ej-i satisfying the conditions of equation ©, respectively. Every term of the computed sequence 
depends on all of the previously computed terms. 

Algorithm 2. Description - Excited State Energies: 

1. Consider A as defined in Algorithm 1. Run Algorithm 1 and let Eo be its output. 

2. Set i = 1 and prepare to compute an estimate of E\. 

3. Repeat steps 1-8 of Algorithm 1, storing the outcome of every measurement. We assume that the 
measurement outcomes are truncated to b bits and we obtain non-negative integers in the range 
{0,..., 2 b — 1}, where b <C t. The role of the extra qubits to — b — t is to increase the success probability 
of QPE. 

4. Take the minimum of the measurement outcomes that exceeds by 2 the last selected outcome and 
mark it selected. This way, with high probability, for each eigenvalue the error will be 0(s ), and the 
algorithm will not produce two different estimates for the same eigenvalue. Note that by taking the 
minimum outcome relative to the previously selected outcome implies that the algorithm does not fail 
to produce estimates for consecutive eigenvalues, unless the eigenvalues differ by O(e). See FigureQ] 

5. Use the values of a and R in the definition of A to rescale and shift the newly selected outcome to 
obtain the estimate E. t . 

6. Set i <— i + 1 and prepare to compute the estimate Ei. 

7. Repeat steps 3 through 6 if i < j. 


01 


02 03 


TO — 1 TO TO + 1 TO + 2 

Figure 1: Example of the selection of the measurement outcomes. Consider three phases fa, fa and fa as 
shown. Assume that the distance between possible outcomes corresponds to error O(e). If to — 1 is selected 
to estimate fa, the next possible outcome the algorithm selects is m+ 1, which provides an estimate for both 
fa and fa in this example. Alternatively, if to is selected to provide an estimate of fa, the next possible 
outcome is m + 2, which provides an estimate of fa, and the algorithm does not care to produce a separate 
estimate for fa. Note that fa and fa are ^-degenerate and either can be ignored. 


8. Output Ei,..., Ej-i. 


It is clear that our procedure as outlined by Algorithms 1 and 2 will produce the j desired estimates JI]) 
of equation (0. However, its cost varies depending on V and the distribution of eigenvalues of L°, which as 
we already mentioned determine the cardinality of S. In the next sections we derive tight estimates for the 
cost and the success probability of our algorithm for particular choices of L° and L. 

We remark that in cases where L and L° are given explicitly, using their properties one may be able 
to obtain a set of trial eigenvectors S with significantly smaller cardinality, substantially improving the 
cost of the algorithm. For example, knowledge of the symmetry groups of L, L°, and V could be used to 
immediately rule out candidate eigenvectors. It is important to observe that different partitionings of the 
Hamiltonian into L° and L — L° may lead to very different sets of trial eigenvectors. Given a Hamiltonian 
L, an important task is to select L° that will result in a relatively small set of trial eigenvectors which can 
be computed efficiently. 


3 Preliminary Analysis 

Consider operators L° and L with the assumptions of Section 12.11 Recall that L° and L are self-adjoint 
operators on a Hilbert Space H with discrete spectra; e.g., see [16]. Then we have the following eigenvalue 
equations: 

L 0 u° = E°u°, i = 0,1, 2..., and let Eg < E? < E§ < ... 

Lui = EiUi , i = 0 , 1 , 2 ..., and let Eq < E\ < E 2 < ■■■ 

Without loss of generality we take all eigenvectors to have unit length. We assume that the eigenpairs 
{Ef , are known. 

First suppose we wish to compute a specific eigenvalue E of L. Let u be a unit vector in the (possibly 
degenerate) subspace associated with E. We have 

II Lu - L°u\\ 2 = || (L° + V)u - L°u\\ 2 = ||Eu|| 2 < \\V\\ 2 L o, 

since u belongs to the intersection of the domains of L° and L. Expanding in the basis of unperturbed 
eigenvectors, we have |u) = JA fa |u°) where fa = (u°|u). Then 

II Lu - L°u II 2 = \\Ej2fa K> - E K> II 2 = II E &(£ - E°) | u°) II 2 

i i i 

Combining these expressions and using the eigenvector orthonormality gives 

||E|| 2 o > ||Lu - L°u \\ 2 = \Pi\ 2 (E°i - E) 2 (10) 
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Assume that there exists c > 1 such that the condition of equation © holds. Observe that this is true for 
instances of the time-independent Schrodinger equation m ■ From equation m, using B > E, we obtain 

Iloilo > £ IAI 2 (E° -E) 2 >J2 m 2 (c\\V\\ L o ) 2 

iei iei 


which we rearrange as 


or equivalently 


D«i 2 A 

iei 

£lAI 2 >i-^=:g 

I 


( 11 ) 


( 12 ) 


Thus there must exist an index k ^ I such that |/3fc| 2 > t^j. If |«S| is not extremely large, then one of the 
first |<S| eigenvectors of L° must have a reasonable overlap)] with u. 


4 Application: Time-Independent Schrodinger Equation 


In this section we consider the time-independent Schrodinger equation on the d-dimensional unit cube with 
Dirichlet boundary conditions to illustrate the algorithms of Section 12.21 that compute the j estimates of 
equation (j4j , where j = 0(1). In particular, consider the eigenvalue problem 

Lu(x) := (— lA + V)u(x) = Eu(x) for all x £ Id ■= (0, l) d , (13) 

u(x) = 0 for all x £ did , 


where V is uniformly bounded by a constant M and has continuous first-order partial derivatives in each 
direction uniformly bounded by a constant C, i.e. || -^rV\\ < C. Thus, without loss of generality we assume 
that V > 0. We set L° = — lA, where 




cP_ 

dxf 


We assume that the eigenvalues of L and L° are indexed in non-decreasing order. We want to approximate 
the first j excited state energies, E ^,..., £’(j_i) , (i.e. the j smallest eigenvalues ignoring multiplicities) 
with error proportional to e, modulo e-degenerate eigenvalues as explained previously. Thus we are interested 
in low order excited state energies because we have assumed that j is a constant. Recall our definitions and 
notation of Section [2] 

As we already indicated, the cost of Algorithms 1 and 2 depends on the cardinality of a set of trial 
eigenvectors S. We will now show that |5| is bounded by a polynomial in d in the case of the Schrodinger 
equation we are considering here. 

The eigenvalues and eigenvectors of L° = — are known to be 


-Eg — -(&i + k 2 + ... + k 2 )n 2 k — (Aq,..., kd) £ 


(14) 


d 

u°j:(x) = 2 d/2 Y\_sm(kinXi) x = (aq,..., x d ) £ [0, l] d k = (Aq,..., k d ) £ N d . 

i=l 

We may re-index them by considering the eigenvalues in non-decreasing order to obtain Eg < E® < ... < 
E® < ... as in Q. Thus Eg = ^dn 2 < i (d + 3)7r 2 = E ° and E ° is a degenerate eigenvalue with dimension 
of its associated degenerate subspace equal to d. Similar considerations apply to the rest of the eigenvalues. 
We remark that the distribution of the eigenvalues of L° is known m- 

We will use <[9|) with c = 2 to derive a set of trial eigenvectors and bound its cardinality. In fact, we derive 
a set of trial eigenvectors that is slightly larger than the set obtained by strictly considering the indices in 

3 Here and elsewhere, by reasonable overlap we mean that the magnitude of the projection is not exponentially small in d. 
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the complement of X in ©. Yet its size is polynomial in d as we will see, and for the sake of brevity, we also 
denote this set by S. In particular, we construct the quantity B of equation ([9]) and show a K = I\(j,V) 
such that for k > K =$■ E® > 2M + B. So we obtain an upper bound for the jth largest eigenvalue of 
L. Clearly the cardinality of S grows with B because we include eigenvectors of L° that correspond to 
increasingly large eigenvalues. The purpose of the construction below is to obtain a crude but helpful in 
our analysis estimate of the distribution of the eigenvalues of L using the eigenvalues of L°; in particular to 
cover possible degeneracy of the eigenvalues of L. 

We select j + 1 values E° s ^ from the strictly increasing sequence of eigenvalues (see ©) such that 


Ko) ■■= E (0) < E (0) +M< E ( S1 ) < E °( S1 ) + M < • ■ • < K,-l) +M< K 


(15) 


where -Eq Sn ) — n = 1) • • • ,3- Indeed it is possible to select a subsequence that satisfies these 

conditions. We know that E° s ^ = \{k\ + k 2 + ■ ■■k d ) 7 r 2 for a certain k. The inequality 




k' 2 + k' m )7T" + ~{k 2 J +1 + fc ^ +2 + ■■■k^)TT 2 > S(° Sn _ 1 ) + M 


is satisfied by selecting m to be a suitable constant and then by selecting k[ > ki + "fi, where 7 i is a suitable 
positive integer constant, i = 1 For example, after fixing to, we can repeatedly increment each 

of the k[, i £ {l,...,m}, successively until the desired inequality holds. Iteratively, we define If?, = 

(k'\ + k'\ + ...k' 2 m )n 2 /2 + (k 2 m+1 + k 2 m+2 + ...k 2 d ) tt 2 /2 for n = 1 , 2 ,.., j. 

By our construction, the interval [.E 1 ^, E® s contains at least j distinct eigenvalues of L , since E^ , — 
E ® s 0 ) > jM, and for every i, E® < Ei < E® + M. Moreover, = E$ + c', where d is a constant. Thus, 
we take c = 2 in ([5]) and define the constant B as 


B:=M + E° (sj) =M + Eq + d. 
From m there exists a K £ N such that 

k>K =» El > 2M + B = 3M + E° (sj) 


(16) 


(17) 


Hence, we construct the set of trial eigenvectors S to be the set of all eigenvectors of L° that correspond to 
eigenvalues less than or equal to 3 M + E® s ,y We bound |S| next. 

The cardinality of S is the number of tuples k £ such that {k 2 + ■ ■ ■ + k 2 ) 7 t 2 /2 < 3 M + dn 2 /2 + d. 
Let to be the number of components k ^,..., ki m of such a k that are greater than 1. Then we have 

(d — to) 7 t 2 /2 + (k 2 ± + ...k 2 m )Tr 2 /2 < 3M + dn 2 /2 + c'. 


Since ki > 2 we have 

3to7t 2 < —mw 2 + (fc 2 + ...fc 2 m )7T 2 < 2(3M + d). 

Hence, to is 0(1). Therefore, in order to construct S one needs to consider tuples k! £ where at most a 
constant number of components are greater than 1. The number of such tuples depends on the number of 
possible combinations by which one can select a constant number of components of k! to be greater than or 
equal to 2. Therefore, this number is polynomial in d@ 

Table Q] below shows the eigenvalues of the Laplacian by considering tuples where a constant number to 
of components exceed 1, assuming that these components are each bounded by a constant N. Observe that 
in all cases, since to = 0 ( 1 ), the multiplicity of the eigenvalues is polynomial in d. 

Therefore, the cardinality of S is polynomial in d. As shown in Section [3l for every eigenvector of L 
that corresponds to an eigenvalue less than or equal to B, there exists an eigenvector of L Q in S such that 
the two eigenvectors have a non-trivial overlap and it follows from (fl2l) that the magnitude squared of this 
projection of the one onto the other will be at least = O(y^y)- 

4 This follows immediately for m = 0(1) from the bound (^) < = poly(d). 
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Table 1: Distribution of eigenvalues of L° = — ^A with respect to the number m of indicies fcj > 2. 


m 

Combinations 

Eigenvalue 

0 


dir 2 / 2 

1 

(?) 

(d — 1)7t 2 /2 + fc l 2 i 7r 2 /2 

2 

(3 

(d-2)ir 2 /2+(kl+kl)ir 2 /2 

l < d 

(?) 

(d-l)r i / 2 + (k* 1 + .. + k*y / 2 


4.1 Finite Difference Discretization 


We obtain a matrix eigenvalue problem by discretizing m on a grid with mesh size h = N € N, using 
finite differences 123112IJj . This yields a matrix Mh := — ^A/j + Vh with size N d x N d . The matrix — ^A/, is 
obtained using a 2d + 1 stencil for the Laplacian {23} p.60]. It is known that the low order eigenvalues of Mh 
approximate the corresponding eigenvalues of L. The eigenvalues and eigenvectors of — h are known and 
are given by 

2 d 

E h,k = ^Esiii 2 (fU,/2) k = (k 1 ,...,k d ) 1 <h<N (18) 

i =1 


U 


0 ^ 
h,k 


where the vectors Vk t £ have coordinates 



(19) 


Vki,i = V2hsin(ki£irh) £ = 1 , 2 ,..., TV % = 1 , 2 ,..., d. 


( 20 ) 


Similarly to (JHJ) , we index the eigenvalues of — h in increasing order ignoring multiplicities to obtain 



E h,( 0) < E h,(l) < ••• 

< E h,(i) < 

(21) 

Then from [39] we have 





K( fc) - E° (k) \ < Cdh 2 

for k = 0(1) 

(22) 

where C > 0 is a constant. 





Vh is an N d x N d diagonal matrix which contains evaluations of V at the grid points truncated to 
|dog 2 /i _1 ] bits of accuracy. Thus Mh is symmetric, positive definite, and sparse. This matrix has been 
extensively studied in the literature [nj uni [23]. For V that has bounded first-order partial derivatives and 
k = 0(1), using the results of [Ml HD] we have that there exists a matrix eigenvalue Eh,k' such that 


| E (k) - E h , k ,\ = 0{dh) (23) 

as dh — > 0, where E^) is defined in ]6|. We will use the algorithms of Section [2721 to approximate the low 
order eigenvalues of Mh , which as we have seen approximate the low order eigenvalues of L. For this, we 
need to construct the set of trial eigenvectors S , and estimate its cardinality. Recall that for the continuous 
operator, the set of candidate eigenvectors is derived using equation and in particular by selecting 

the eigenvectors of L° that correspond to eigenvalues less or equal to 2 M + B = 3 M + E® s . ^. So for the 
discretized case we select the eigenvectors of — |A h that correspond to eigenvalues less than or equal to 
3 M + E? , + 0(dh 2 ) due to equation (1221) . Since dh —> 0, without loss of generality we slightly modify 
equation El, to select the eigenvectors of Mh that correspond to eigenvalues less than or equal to 

2M + B = 3M + E ° {s . } + 1 (24) 

for sufficiently small h, where this equation effectively redefines B by increasing its value by 1. Thus the 
cardinality of S in the case of the matrix Mh follows from the continuous case and remains polynomial in d. 
Specifically, we define 

5 := K lfe : E° h ,k < 3M + E° (sj) + 1} (25) 
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4.2 Algorithm for Excited State Energies 

We now give the details of Algorithms 1 & 2 of Section 12.21 applied to the time-independent Sclirodinger 
equation m- Given e, the algorithms produce the j eigenvalue estimates E 0 < ... < Ej_ 1 of equation ©• 
Algorithm 1 computes Eq■ For this, QPE [26] is applied repeatedly with its initial state taken to be every 
single element of the set of trial eigenvectors S. We use repetitions of the procedure to boost the success 
probability. We remark that our Algorithm 1 computes the ground state energy in a way similar to 129] [28] , 
but under weakened assumptions. Algorithm 2 iterates j — 1 times the procedure of Algorithm 1, at each 
iteration producing the next estimate Ei by taking into account all the previously produced estimates as we 
will explain below. 

Both algorithms use QPE as the main module. The purpose is to compute approximations of the 

eigenvalues of the matrix Mh of the previous section. Setting N = 2 T 2 lo S 2 _ we discretize m with mesh 

size h = -jyFj < to obtain the N d x N d matrix Mh, where we have N d = 0((j) 2d ). From C51) . we obtain 

that the low order matrix eigenvalues approximate the low order eigenvalues of the continuous operator with 

2 

error proportional to dh = 0(^-). The reason we have taken very small h is because we want to ensure that 
e-degenerate eigenvalues of the continuous operator will be approximated by tightly clustered eigenvalues 
of Mh- As M is a constant, without loss of generality we may assume that e -1 M. Since the largest 
eigenvalue of —h is bounded from above by 2 dh~ 2 , and V is uniformly bounded by M , we obtain that 
\\M h \\ is bounded from above by 2 dh~ 2 + M <C 3 dh~ 2 , in the sense that = o(l). 

Let R = 3dh~ 2 and consider the matrix W = e zMh / R . Its eigenvalues are e lEh ^ R = e = e 2m ^, where 
Eh is an eigenvalue of Mh and (f> := 4=^ denotes the corresponding phase. 

QPE is used to compute an approximation 0 of 0 with b = 5[log 2 yl + 7 bits of accuracy, and from this 
we get E = 2'KR<j) so that 

\E-E\<\E- E h \ + \E h -E\ = 0(e) (26) 

where E denotes the eigenvalue of L that Eh approximates according to (fl3l) . QPE uses two registers, the 
top and the bottom. The size of the top register is related to the accuracy of QPE and its success probability. 
Recall that QPE succeeds when it produces an estimate with accuracy 2 -b . The bottom register is used 
to hold an (approximate) eigenvector of Mh corresponding to the phase of interest, and therefore has size 
dlog 2 N = d - 0(log d ). The number of qubits in the top register is t = b + t 0 , so that QPE has accuracy 2~ b 
with probability at least 1 — 9 ( 2 *o_ 2 ) > assuming that an exact eigenvector is provided as initial state in the 

bottom register [221 Sec. 5.2]. QPE uses powers of W, namely W 2 °, W 2± ,..., W 2± ~ 1 . We will approximate 
these powers using a splitting formula with error, as we will see below. This reduces the success probability 
of QPE to at least p := 1 — 2<0 1 _ 2 . We will set to to be logarithmic in d, and will give all the details later on 
when dealing with the cost of our algorithm. 

Consider an eigenvalue Eh < B + 2M (see equations © and (l24l) l of the matrix Mh and let Uh denote 
an eigenvector corresponding to Eh- Then QPE with initial state some i £ S succeeds with probability 
at least p Uh (i) := !«,/ -p = • (1 - ^ 2 ) DQ- 

Recall that S contains eigenvectors of L° that correspond to eigenvalues E® < B + 2M as defined in (12511 , 
and that the cardinality of S is polynomial in d. Applying the same approach of Section[3]for the eigenvectors 
of Mh, we conclude that for every eigenvector Uh of the matrix Mh that corresponds to an eigenvalue less 
than B, there exists a vector k £ S such that \u R u° hk \ 2 > where we have used equation m with 
c = 2 (since the value of B we are using here leads to c = 2 in this case too). Thus, after we run QPE with 
each element of S as initial state, the probability that at least one of the outcomes (in principle we do not 
know which one) will give a good estimate of Eh is at least p Uh (k) > 

We repeat the whole procedure r times to boost the success probability of obtaining an estimate of E h 
with accuracy e. Indeed, the probability that QPE fails with all initial states taken from S and in all its 
r|«S| repetitions is 

|li(l-Pu,C0)j <(1 -Pu h (k)) r <e- r ^<e- r *fc p (27) 

Thus, the probability that at least one of the r|<S| outcomes will lead to an approximation of Eh with accuracy 
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0(e) is at least 


1 _ e - r 4]%T p = 1 _ e r WT' (1 To 33 ) (28) 

We can boost this probability to be arbitrarily close to 1 by taking r = poly(d), since |<S| is polynomial in d. 

Observe that Algorithm 1 selects the minimum measurement outcome from all the runs of QPE, and uses 
it to obtain Eq. Let this outcome be m' G {0,..., 2* — 1}. The algorithm converts m! to mo = \_m'2~ to \ G 
{0,..., 2 b — 1} and uses it to obtain E 0: according to the formula 

E 0 = 2TrRj> 0 = 2t tR^-. (29) 


Since |m , /2‘ — 4>q\ < 2~ b and (j)o,m' /2* G [mo/2 b , (mo + l)/2 fa ], it follows that |27ri?<(>o — Eq\ = “2t*R • \<fio — 
7 Tio/ 2 b | < 2irR/2 b < e, which together with (E51) gives (E()l) . 

Let Go = {m : |p- — (f> o| < 775 -}. Algorithm 1 fails either if none of the converted outcomes is an element 
of Go, or at least one of the converted outcomes is an element of Go, but there is another converted outcome 
(produced by a failure of QPE) smaller than the minimum element of Go- Thus, we can bound the total 
probability of failure by 


Pr(Algorithm 1 fails) 


+ 

< 

+ 


< 


< 

< 

< 


Pr(none of the outcomes leads to an element of Go) 
Pr(one of outcomes leads to an element of Go 
but there is at least one other smaller converted 
outcome) 

e -r 4j%rP 

Pr(QPE failed in at least one of 

the rj«S| runs) 

e~ r ^m p 

(1 — Pr(every run of QPE approximates 
one of the phase with error 2 ” 6 )) 

e ~ r im p + (1 - j/l 5 l) 


—r 

e 


_ 2 _n- 
41 «s | ^ 



e r 4|s |( 1 


2 4 ° - 2 ’ 


(30) 


where the third from last inequality follows from equation m below. Observe that this bound can be made 
arbitrarily close to 0 by selecting the number of repetitions r to be a suitable polynomial in d, since |<S| is 
polynomial in d, and by taking to = /31og d, where /3 is an appropriately chosen constant. We have used the 
fact that if a measurement outcome £ fails to estimate any of the phases, i.e. — (f> s \ > 2~ b for all phases 
<f> s corresponding to eigenvalues of M f,, then 


N —1 


Pr (^) = 22 c sl a (^’^)l 2 - 

s=0 


N —1 

E 

s=0 


1 


1 


2*o _2 2 * 0-2 


= (1 — p)- 


(31) 


Here, the c s denote the projections of the initial state onto each of the eigenvectors of M/,., and the \a(£, 4> s )\ 2 
denote the probability to get outcome £ given the exact eigenvector Uh lS as input. We have upper bounds 
for these quantities from [2B1 Eq. 5.34]. Therefore, the probability the measurement outcome estimates at 
least one (or, some) phase is 1 — Pr{£) > p. 

Recall that the set S has been constructed using an upper bound for E(j_iy, see equations (flSl) and 
m- Algorithm 2 essentially repeats Algorithm 1 (j — 1) times, but selects the converted measurement 
outcome in a different way by considering the already selected outcomes. At repetition i. it selects the 
minimum converted outcome mi that exceeds the outcome selected at the previous iteration by at least 2 , 
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i.e. rrii > rrii-i + 2, where rrii = |_to( 2~ to \ and to' is a measurement outcome at the zth run, i = 1,2, 3, j — 1; 
see also equation (1291) . The success probability for both Algorithm 1 and Algorithm 2 follows from (l30l) and 
is at least 


1 - 




(32) 


which can be made arbitrarily close to 1 by selecting r to be a suitable polynomial in d and taking to to be 
sufficiently large. 

Note that Algorithm 2 computes Ei = 2i i = 1 , 2 ,..., j — 1 , as estimates of the eigenvalues according 
to equation (0) and the conditions Cl and C2 that follow it. If both algorithms are successful with high 
probability, at the zth run we have that there exists a phase 4>i corresponding to an eigenvalue of Mh such 
that |27 iR(j)i — Ei\ < < e. The condition nii > m.j_i + 2 in the selection of measurement outcomes 

guarantees that for any two i\ ^ Z 2 , the computed matrix eigenvalue approximations satisfy E q ^ Ei 2 and 
| E^ — Ei 2 1 = fl(e) because for the corresponding phases we have (j)^ ^ cf>i 2 as belonging to different intervals; 
see Figure HJ Moreover, the E ix and E i2 also approximate different eigenvalues E ix ^ E i2 of the continuous 
operator because we have used a very fine discretization. Finally, the algorithm does not fail to produce 
consecutive eigenvalues unless they differ by less than 0(e) because we always select the minimum outcome 
that satisfies rrii > rrii -1 + 2. 


4.2.1 Cost of Quantum Phase Estimation 

Algorithms 1 & 2 use QPE as a module. The cost of QPE depends on the cost to prepare its initial state, and 
on the cost to implement the matrix exponentials W 2 ,W 2 , ..W 2 , where W = e lMh ^ R . We approximate 

these exponentials below using Suzuki-Trotter splitting, the analysis of which proceeds similarly to that of 

E3IH 

The initial states are taken from S which contains eigenvectors of — \Ah according to (l?5l) . Each eigenvec¬ 
tor can be prepared efficiently using the quantum Fourier transform, which diagonalizes the Laplacian, with 
a number of quantum operations proportional to d-log 2 ^ and using number of qubits log 2 N d = d- 0(log ^). 
We remark that from the tensor product structure of the eigenvectors of — f A h, it suffices to prepare eigen¬ 
vectors of the one-dimensional Laplacian; see e.g. 011111113- 

Now let us turn to the approximation of the matrix exponentials. We simulate the evolution of the 
Hamiltonian H = Mh/R for times 2 r , r = 0,1,..., t — 1, where we have set t = b + to- Let H = H i + H -2 
where Hi = — A/,/21?, and H 2 = Vh/R, where we assume V is given by an oracle. 

To simulate quantum evolution by Hi, assuming the known eigenvalues of — ^A^ are given by a quantum 
query oracle with 0(log i) bits of accuracy, we again use the quantum Fourier transform to diagonalize Hi 
with cost (i.e., a number of quantum operations) bounded by d- 0( log 2 d ), and requiring a number of qubits 
proportional to dlog 7 . Alternatively, if the eigenvalues of -|A(,, are implemented explicitly (without an 
oracle) by the quantum algorithm, then the number of quantum operations required is a low order polynomial 
in d and log 2 -, and so is the number of qubits 0. For simplicity, we will not pursue this alternative here. 
The evolution of a system with Hamiltonian H2 can be implemented using two quantum queries returning 
the values of V at the grid points, and phase kickback. The queries are similar to those in Grover’s algorithm 
j26l and the function evaluations of V are truncated to 0(log y) bits. 

We use a splitting formula S 2 k of order 2k + 1, k > 1, to approximate W 2 = e l ( H i+ H 2 ) 2 by a product 
of the form 

N t 

J\e iAtZl , (33) 

1= 1 

where An £ {Hi, # 2 } and suitable Z( that depends on t and k. 

The splitting formula S- 2 k is due to Suzuki [351136] . It is used to approximate e d B + c ) At ; where B and C 
are Hermitian matrices. This formula is defined recursively by 

S 2 {B,C,At) = e iBAt/ 2 e iCAt e iBAt /2 

S 2k (B,C,At) = [S2 k -2(B,C,p k At)] 2 S2k-2(B,C,(l-4p k )At) 

x [S2k-2{B, C , PfcAf)] 2 , 
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where pk = (4 — 4 1 /( 2fc - 1 )) 1 , fc = 2,3,.... 

Unfolding the recurrence above and combining it with m Thm. 1] we obtain that the approximation of 
W 2 has the form 

pj/2 T _ e iHia.T t o e iH2b-r t i e iHia r ,i . . . e iH 2 b Tt L T e iH 1 a r ,L T (34) 


where a Tj o, • • •, a r ,L T and 6 T> i, ..., & r ,L r and L T are parameters, r = 0, ... ,to + b — 1. The number of 
exponentials involving Hi and H 2 in the expression above is N r = 2 L T + 1. An explicit algorithm for 
computing each W 2 is given in | |29] . 

Let || • || be the matrix norm induced by the Euclidean vector norm. From pTTl Thm. 1 & Cor. 1] the 
number N t of exponentials needed to approximate W 2 by a splitting formula of order 2fc + 1 with error e T , 


t = 0,... ,to + b — 1, is 


N T < 16e||Fi||2 r 


py- 1 ^ 8e 2 T || JL 2 II 


1/(2 fe ) 


for any k > 1. Since we want to approximate all the W 2 , r = 0,1, ...,to + b — 1, we sum the number 
of exponentials required to approximate each one of them. Thus the total number of matrix exponentials 
required by Algorithm 2, Mot, is bounded from above by 


to~\-b—l 

N tot = jr\S\ N r 

T —0 

/ /oc\ k—1 tp -\-b— l / or \ 1/(2A:)' S 

< jr\S\ l 16e||M|| (8e||tf 2 ||) 1/(2fc) £ T 


(35) 


product r\S\/g(d) = o(l) (as d —► 00 ). We then select the error of each exponential to be e r = 
r = 0,... ,io + b — 1. It is easy to check that 


jt—O — 20g(d) ‘ 


The factor jr|5| is the number of executions of QPE performed by our algorithms, and the second factor is 
the cost of a single QPE. Note that j is the number of eigenvalues we wish to estimate, |<S| is the number 
of eigenvectors we use as initial states, and r is the number of times we repeat QPE per initial state to 
boost the success probability of getting the desired outcome. We select a polynomial g{d) such that the 

2T+1 —(b+to) 
40 W) ’ 

Thus the success probability of 
QPE is reduced by at most twice this amount [26] P- 195], giving 1 — 2 ( 2 *o- 2 ) — lOg(d) ■ Nexb we set 
to = |k>g 2 ( 55 (d) + 2)J, to get p= 1 — 2 tg_ , that we used above in deriving equation (I3U1) . Our choice of g(d) 
and to aims to make the bound of equation (1321) arbitrarily close to 1. 

The largest eigenvalue of —A h is Adh~ 2 sm 2 (nNh/2) < Adh~ 2 . Since R = 3 dh~ 2 , Hi = — ^A h/R = 

— I gdF 2 ^ h an< ^ we li ave 11 Mil < = §• Since V is uniformly bounded by M and H 2 = Vh/Rwe 

have |jif 2 || < M/3dh~ 2 . Substituting the value of e T in (1551) . yields that the algorithm uses a number of 
exponentials of Hi and H 2 that satisfies 


Mot < jr\S\ ( 16e||M|| (8e||Id2||) 1/(2fc) J £ 2 " ( 

< jr\S\ (l6e\\Hi\\ (8e|| ||) 1/(2fc) j (20 ff (d)2‘°+ b ) 
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Using the bounds on ||iU|| and \\H 2 W, we obtain 
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From b = 5[log 2 f ] + 7, we have 2 b = 2 5 ^ll+ 7 < 2 12 (f) 5 = 0(|(). Since h< we have 2 b /i 2 < 2 12 f = 
0(f). Also, 2 to = 2L 1 o S2(5sO)+2)J < 5 5 (d) + 2 . We obtain 


M tot < j>|5| (^(5 5 (d) + 2)(2 1 


x ( 160e (5 ff (d)+2) ( 2 12 ^j 


d 

£ 

d\ M 


k- 1' 


l/(2fc) 


< J>|5| O 


d 5 ff (d) /^25y ^gi c/ 2 (d)^ 


l/(2fc)N 


(36) 


for anv 
The 
is given 


k > 0, where C and C are suitable constants. 

optimal k*, i.e., the one minimizing the upper bound for fiftot in (ED: 

by 


k* 


9 l°§25/3 



1 

2 


= C\ In 


is obtained in El Sec. 5] and 


for a suitable constant C, since g(d) is a polynomial in d and we are taking its logarithm. With k* and using 
again El Sec. 5], equation ED yields 



O g 2 (d) 



as de —> 0, 


(37) 


where we have used j = 0(1) and r|<S| = o(g(d)), and where the equality above holds asymptotically for 
arbitrarily small rj > 0. 

We remark that of the JV* ot matrix exponentials roughly half involve H\ and the remaining involve H 2 ] 
see ED- Since each exponential involving H 2 requires two queries the total number of queries is also of 
order N^ ot . The cost to prepare the initial state, to diagonalize — |A/j, and to implement the inverse Fourier 
transform that is applied prior to measurement in QPE, is proportional to 

dlog 2 ^ + (t 0 + b ) 2 = O ^d log 2 y^ , 


since to + b = 0(log -). Hence, the total number of quantum operations, excluding queries, is proportional 
to 

K*-dl Og 2 -. (38) 

£ 

Equations (1371) and (1551) yield that the total cost of the algorithm, including the number of queries and the 
number of all other quantum operations, is proportional to 


dg 2 (d) 



where S > 0 is arbitrarily small. Finally, using equation (1321) we can select r to be polynomial in d and 
obtain success probability at least |, and the cost remains polynomial in - and d. We summarize our results 
in the following theorem. 


Theorem 1 Consider the time-independent Schrodinger equation m on the d-dimensional unit cube with 
Dirichlet boundary conditions and where the potential V and its first-order derivatives are uniformly bounded. 
Algorithms 1 & 2 of Section \2.2\ comvute approximations of j = 0(1) low order eigenvalues with error 0(e), 
and satisfying conditions Cl and C2 of Section^ with probability at least 



2*o - 2 


17 












where r and |<S| are polynomial in d, to = Llog 2 (5g(d) + 2)J, and g(d ) is a polynomial in d selected such that 
r|<S| = o(g(d)). The algorithms apply QPE with initial state each element of a set of trial eigenvectors S , 
and repeat this procedure r times. 

They use a number of queries proportional to 

5+5 

g 2 (d) as de -+ 0 , 
a number of quantum operations excluding queries proportional to 

5+5 

d g 2 (d ) as de -+ 0 , 

where 5 > 0 is arbitrarily small. The algorithms use a number of qubits proportional to 

d log - + log g(d). 

£ 




Remark 1 The 5 in the exponent of d is due to the fact that for simplicity we have taken N = 2 ^ log2 i 1 in 
the discretization of the continuous operator and our consequent choice of b, the number of bits of accuracy 
of QPE. As we explained, the purpose of the fine discretization is to ensure that degenerate eigenvalues of the 
continuous problem are approximated by tightly clustered eigenvalues of the matrix. It is possible to reduce 
the cost estimates of Theorem Q] by taking N = log2 si and b = (3 + 27 ) [~log 2 + 7 where 7 G (0,1) is 

a suitable constant. Then the total cost becomes proportional to d g 2 (d) (^) 3+<5+ “ 7 p We d 0 no i pursue this 
any further since out goal was to establish an algorithm with cost polynomial in d and £ _1 . 

On the other hand, the classical complexity of approximating a constant number of low order eigenvalues 
with error e grows as (i) in the deterministic worst case. Therefore, the problem suffers from the curse of 
dimensionality. The lower bound follows from the corresponding lower bounds for approximating the ground 
state energy with ||V|| < 1 , since more general conditions on V are considered in this paper. A detailed 
discussion concerning the classical complexity lower bounds for approximating the ground state energy can be 
found in flf. Since our quantum algorithm for this problem has cost polynomial in d and it vanquishes 
the curse of dimensionality. 


5 Conclusion 

There are a number of recent results suggesting that certain eigenvalue problems are very hard, even for 
quantum computers. On the other hand, obtaining positive results for eigenvalue problems showing advan¬ 
tages of quantum computers over classical computers is particularly important. We discuss such a positive 
result for the approximation of the ground state energy in 1301 . yet there is much more to be done. 

In this paper, we consider the approximation of ground and excited state energies (low order eigenvalues) 
of a self-adjoint operator L. Typically, L is discretized to yield a matrix eigenvalue problem. Since L is self- 
adjoint, the resulting matrix is symmetric. It is important that the discretization is such that the eigenvalues 
of interest of L are approximated accurately by matrix eigenvalues. This usually increases the matrix size. 
There are numerous classical algorithms that approximate matrix eigenvalues and/or the corresponding 
eigenvectors. In the case of symmetric matrices, we can approximate the eigenvalues that belong to a given 
range using the bisection method. In general, the cost of classical algorithms is bounded from below by the 
matrix size. Thus, the problem becomes hard when the matrix size is huge. On the other hand, quantum 
algorithms provide a way of overcoming this difficulty in certain cases. This can be accomplished using 
QPE as a module for approximating individual eigenvalues. QPE is efficient as long as two conditions are 
met. The first condition is that we are able to perform efficient quantum Hamiltonian simulation for the 
matrix whose eigenvalues are sought. There are numerous papers in the literature providing conditions 
and algorithms for efficient Hamiltonian simulation on a quantum computer. When a Hamiltonian is given 
by an oracle, efficient simulation can result from quantum parallelism and splitting formulas. The second 
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condition needed in QPE is that we can prepare efficiently a quantum state encoding an approximation of 
an eigenvector corresponding to an eigenvalue of interest. It is enough that the approximate eigenvector 
is relatively good but by no means perfect. It suffices that the approximate eigenvector does not have an 
exponentially small (in the problem parameters) overlap with the unknown eigenvector. In our case, since 
we will be approximating a number of eigenvalues, we need an equal number of approximate eigenvectors. 
Obtaining the necessary approximate eigenvectors is a difficult task for arbitrary L. 

We propose a way to overcome this difficulty and construct a set of approximate eigenvectors, which we 
call trial eigenvectors , using the structure of L. Since L = L° + V, and we have assumed that we know the 
eigenvalues and eigenvectors of L , we use a perturbation argument to construct the set of trial eigenvectors. 
We select eigenvectors of L° corresponding to a particular range of its eigenvalues. Equivalently, we exclude 
the eigenvectors of L° that correspond to eigenvalues which significantly exceed the eigenvalues of L that we 
wish to approximate. We describe an algorithm that uses QPE, and the set of trial of eigenvectors, in order 
to approximate low order eigenvalues. If the cardinality of the set of trial eigenvectors is relatively small, i.e. 
its size is at most polynomial in the problem parameters, then our algorithm is efficient. 

In summary, general conditions for the approximation of ground and excited state energies on a quantum 
computer follow by combining conditions for efficient quantum simulation and for deriving a relatively small 
set of trial eigenvectors that can be implemented efficiently as quantum states. 

We illustrate how these general conditions are met in a special case of the time-independent Sclirodinger 
equation with d degrees of freedom. We show how our algorithm approximates a number of low order 
eigenvalues with error e and high probability. From our earlier work on this problem, we know that the 
complexity of approximating the ground state eigenvalue on a classical computer grows exponentially with 
the number of degrees of freedom in the worst-case. Therefore, under the same or weaker assumptions as in 
this paper, the approximation of low order eigenvalues on a classical computer satisfies the same lower bound. 
The problem suffers the curse of dimensionality in the classical worst case. For quantum algorithms, our 
previous approaches for computing the ground state energy require stronger conditions on V than those we 
consider here, and these approaches do not extend to computing excited state energies. We have developed 
an entirely new approach to approximate not only the ground state energy, but also excited state energies, 
with cost polynomial in d and e . Our quantum algorithm vanquishes the curse of dimensionality. 

We remark on several open problems. We have assumed that L = L° + V, but such a partition need not be 
unique, and different partitions may result in algorithms with significantly different costs. It is possible, that 
with additional assumptions, one would be able to determine suitable partitions leading to fast algorithms. 
Such a characterization is an open problem. We have provided a condition for constructing a set of trial 
eigenvectors S. Improving this condition to minimize the size of the resulting set S is another open problem. 
Finally, in our initial investigation of quantum algorithms for eigenvalues problems, we considered strong 
assumptions on V and obtained efficient algorithms. Progressively, we have weakened these assumptions. It 
is important to continue working in this direction to extend the scope of our algorithm. 
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